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Abstract:  The  problem  is  to  calculate  an  approximate  solution  of  an  initial  value 

problem  for  a scalar  autonomous  differential  equation.  A generalized  notion  of  a 
nonlinear  Runge-Kutta  (NRK)  method  is  defined.  We  show  that  the  order  of  any 
s-stage  NRK  method  cannot  exceed  2s  - 1;  hence,  the  family  of  NRK  methods  due  to 
Brent  has  the  maximal  order  possible.  Using  this  result,  we  derive  complexity  bounds 
on  the  problem  of  finding  an  approximate  solution  with  error  not  exceeding  s.  We  also 
compute  the  order  which  minimizes  these  bounds,  and  show  that  this  optimal  order 
increases  as  • decreases,  tending  to  infinity  as  i tends  to  zero.  . 
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1.  Introduction 


Let  2>  be  a subset  of  the  real  numbers  R,  and  let  V ■ {v  : domain(v)  cR-*RJ 
be  a set  of  functions,  such  that  the  initial  value  problem  of  finding  a function 


x : [0,  1]  R satisfying 


«(t)  - v(x(t))  0 < t < 1 


x(0)  - Xq 

has  a unique  solution  for  every  (xq,v)«  5)xty  (The  differential  equation  in  (1.1)  is 
said  to  be  a scalar  autonomous  differential  equation.)  We  are  interested  in  the 
computational  complexity  of  using  one -step  methods  to  generate  an  approximation  to 
(1.1)  on  an  equidistant  grid  (in  the  sense  of  Stetter  [73]h  that  is,  the  methods 
considered  give  approximations  Xj  to  x(ih)  by  the  recurrence 

(1.2)  xj+1  - Xj  + h *(Xj,h)  (0  i i S n - 1), 

where  h - n-*  is  the  step-size  of  a grid  with  n points,  and  f is  the  increment  function 
(Henrici  [62])  for  the  method.  (For  brevity,  we  will  refer  to  "the  method  f") 

In  Werschulz  [76a],  we  discussed  the  complexity  of  solving  autonomous  pvstems 
of  differential  equations;  in  this  paper,  we  will  consider  only  the  case  of  a single  scalar 
autonomous  equation.  Clearly,  the  results  of  Werschulz  [76a]  hold  for  problems  of  the 
form  (1.1).  However,  in  this  paper  we  will  discuss  the  complexity  of  solving  (1.1)  via 
nonlinear  Runee-Kutta  methods  (abbreviated,  "NRK  methods").  We  only  consider  the 
scalar  case  (1.1),  since  it  is  not  known  whether  NRK  methods  exist  for  more  general 
systems. 

In  Section  2,  we  give  the  formal  definition  of  "NRK  method,"  and  show  that  no 
NRK  method  using  s evaluations  of  v ("stages")  can  have  ordpr  exceeding  2s  - 1.  Thus, 


the  set  of  s-stage  methods  of  order  2s  - 1 described  in  Brent  [74]  has  maximal  order 
in  the  class  of  NRK  methods. 

In  Section  3,  we  use  the  results  of  Brent  [74]  and  Section  2 to  find  upper  and 
lower  bounds  on  the  complexity  of  finding  an  approximate  solution  whose  error  does 
not  exceed  «,  using  a method  of  fixed  order.  These  results  are  then  used  to  calculate 
optimal  orders  which  minimize  these  complexity  bounds.  We  show  that  the  optimal 
order  increases  as  * decreases,  tending  to  infinity  as  • tends  to  zero.  Finally,  we 
compare  the  complexities  of  NRK  methods,  Taylor  series  methods,  and  linear  Runge- 
Kutta  methods.  We  show  that  the  best  NRK  methods  Known  are  asymptotically  better 
(as  • tends  to  zero)  than  the  best  linear  Runge-Kutta  methods  possible,  but  are 
asymptotically  worse  than  the  best  Taylor  series  methods  Known  If  the  cost  of 
evaluating  the  derivative  of  v is  bounded  for  all  K. 


2.  Maximal  Order  for  NRK  Methods 


Before  proceeding  any  further,  we  will  review  some  basic  notions  from 
Werschuiz  [76b}  The  following  notational  conventions  will  be  used.  Let  X be  an 
ordered  ring;  then  "X+"  and  "X'M'"  respectively  denote  the  nonnegative  and  positive 
elements  of  X.  (This  is  used  in  the  cases  X - R,  the  real  numbers,  and  X ■ Z,  the 
integers.)  The  symbol  means  "is  defined  to  be."  We  use  "I"  to  denote  the  unit 
interval  [0,  1}  The  notations  "x  A a"  and  "x  t a"  are  used  to  indicate  one-sided  limits, 
as  in  Buck  [65}  Finally,  if  xj  * X2  : R **  and  ea : -♦  R are  differentiable,  then 

for  i - 1,  2,  we  write 

dj  »(xi(t),  X2^) 

for  the  result  of  differentiating  «(xj,  X2^  w'^  r®sPect  to  Xi>  and  then  substituting 
Xi  ■ Xi<t),  X2  " X2<t)- 

We  next  describe  the  model  of  computation  to  be  used.  We  assume  only  that  all 
arithmetic  operations  are  performed  exactly  in  R (i.e.,  infinite-precision  arithmetic)  and 
that  for  all  v < 19,  we  are  able  to  compute  the  value  of  v at  any  point  in  its  domain.  In 
addition,  we  must  pick  an  error  measure,  so  that  we  may  measure  the  discrepancy 
between  the  approximate  solution  produced  by  p (via  (1.2))  and  the  true  solution.  For 
the  sake  of  definiteness,  we  use  the  global  error 

(2.1)  max  OsiSn  lx(ih)  ' xil  • 

Other  error  measures  may  be  used,  such  as  the  local  error  BSL  step  and  the 
error  per  unit  step  (see  Henrici  [62]  and  Stetter  [73]  for  definitions);  this  would 
involve  only  a slight  modification  of  the  results  contained  in  the  sequel. 

Finally,  we  will  say  that  * - {*p  : p « Z 4"f}  is  a basic  sequence  of  methods  if 
there  exist  functions  « : R+xl  -»  R and  , «u  : R + •*  R*  such  that 
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(2.2)  *G(*p'h>  “ "<p,h)  hP  (or  h i I and  p € Z +*  , 

where 

(2.3)  0 < «L(p)  S «(p,h)  S «u(p)  < +«  for  h « I . 

We  say  that  *p  has  order  p.  This  is  a slight  extension  of  the  definition  of  order  given 
in  Cooper  and  Verner  [721  the  function  introduced  here  is  necessary  and  sufficient 
for  the  "order"  of  a method  to  be  unique.  (Here  we  introduce  the  convention  of 
attaching  the  subscripts  "L"  and  "U"  to  quantities  dealing  with  lower  and  upper  bounds 
(respectively)  on  complexity.) 

We  now  consider  a generalization  of  the  familiar  linear  Runge-Kutta  methods 
which  are  found  in  standard  texts  such  as  Henrici  [62].  A basic  sequence  • is  said  to 
be  a sequence  of  nonlinear  Runge-Kutta  methods  ("NRK  methods")  if  each  increment 
function  *p  « • m*y  be  written  in  the  form 

(2.4)  *p(Xj,h)  rs(xQth;  K0, ...  .Kg.p  , 

where 

(2.5)  kj  v(yj) , yj  fj(Xj,hi  k0, ...  .kj.j)  (0  S j S s - 1) 

for  suitable  functions  tj  : RxRxrI  ■*  R (0  S j £ s).  We  say  that  ^p  has  s • s(p) 
stages,  so  that  an  s-$tage  NRK  method  uses  s evaluations  of  v.  Since  the  one-step 
method  ^p  defined  by  (2.4)  and  (2.5)  is  stationary  (i.e.,  does  not  change  from  step  to 
step),  we  need  only  describe  how  xj  is  generated  from  xq. 

Brent  [741  [76]  considered  the  problem  of  finding  a simple  root  f of  a nonlinear 
function  FsR  -*F,  using  the  Brent-information  (Meersman  [76]) 

(2.6)  Wb#8(F)  {F(x0),F'(x0),F'(y1),...,F'(ys_1)}, 

where  xq  is  an  initial  approximation  to  f,  and  yj  , ... , ys_i  are  to  be  determined.  Let 
x j be  a sufficiently  good  approximation  of  the  appropriate  zero  of  the  minimal-degree 
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polynomial  interpolating  the  information  92giS(F).  Then  Brent  [74]  showed  how  to 
choose  , y8_j  so  that 

(2.7)  l*i  -fl  - CXb.0  -f|2s>  asx0-r. 

This  defines  an  iterative  method  of  order  2s  for  finding  f. 

Let  us  now  define  a function  F by  setting 

(2.8)  F(z)  f2  d|/v(|)  - h , 

*0 

and  note  that  x(h)  is  the  zero  of  F.  Recalling  that  order  for  iterations  is  defined 
differently  than  is  order  for  one-step  methods,  (2.8)  shows  how  an  s-stage  NRK 
method  of  order  p may  be  derived  from  a (p  ♦ l)*h-order  iterative  method  for  zero- 
finding which  uses  the  Brent-information  (2.6).  Using  this  transformation  and  (2.7), 
Brent  [74],  [76]  exhibited  a sequence  #mqRK  of  "m odified " Brent-Runge-Kutta  methods 
(MBRK  methods"),  in  which  the  s-stage  method  has  order 

(2.9)  p - 2s  - 1 . 

Furthermore,  Meersman  [76]  proved  that  this  order  is  the  greatest  possible  in  the 
class  of  all  such  BRK  methods.  We  now  extend  Meersman’s  result  to  include  all  NRK 
methods. 

Theorem  2.1:  No  s-stage  NRK  method  can  have  order  greater  than  2s  - 1. 

Proof:  Let  f be  an  s-stage  method  with  order  p.  We  will  construct  (from  p)  an 
iterative  method  if  of  order  q :•  p ♦ 1 for  finding  a simple  zero  f of  an  arbitrary 
analytic  function  F : R -*  R. 

The  method  + is  defined  as  follows.  Let  xq  be  an  approximation  to  f such  that 
F'  is  nonzero  between  xq  and  f.  (Since  F'(f)  t 0,  such  an  xq  exists.)  Write  (q  :•  F(xq)i 
without  loss  of  generality,  assume  tQ  < 0.  Now  apply  one  step  of  ?,  using  a step-size 
of  -tQ,  to  the  problem 

*(t)  - F'<x(t))"1  <t0  < t < 0)  with  x(t0)  - x0, 
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(whose  solution  is  the  functional  inverse  of  F,  so  that  x(O)  - F l(0)  • f)»  then  # Is 
given  by 

*(x0)  !"  x0  “ *0  ^x0'_t0*  • 

By  the  definition  of  order  for  iterative  methods,  it  is  clear  that  # has  order  qi 
moreover,  * uses  the  generalized  Brent  information  (Definition  IL3.8  of  Meersman  [76]) 

Wgb,s(F)  {F(xo>.F,(yo>,F'(y1),...,F/(ys.i»  • 

Suppose  that  yo  + x0;  then  q s 2s  by  Theorem  II.3.3  of  Meersman  [76}  On  the  other 
hand,  if  yo  - x0,  then  i uses  the  Brent-information  (2.6);  by  Theorem  IL2.4  of 
Meersman  [76]  (also  due  to  Woiniakowski),  we  have  q S 2s  in  this  case  also.  Thus  in 
either  case,  we  find  that 

p ♦ 1 • q S 2s  , 
and  the  desired  result  follows.  | 

Thus  is  intormationallv-optimal  in  the  class  of  NRK  methods,  in  the  sense 

that  each  in  *mBRK  uses  the  ni'nimum  number  of  stages  possible  for  a p^-order 
NRK  method. 
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3.  Complexity  Bounds  for  NRK  Methods 

In  this  Section,  we  will  compute  lower  and  upper  bounds  on  the  total  number  of 
arithmetic  operations  C(p,a)  required  to  guarantee  that  if  is  a p^-order  NRK 
method,  then 

(3.1)  frgVp.h)  £ i e'* 

for  a given  p f I**  and  a ( H++.  (Here  e is  the  base  of  the  natural  logarithms.)  Since 
a > 0,  we  have  0 < < < 1;  clearly  a increases  as  i decreases,  and  « tends  to  infinity  as 
• tends  to  zero. 

In  the  methods  we  consider,  we  may  write 

(3.2)  C(p,«)  - n c(p)  - h"1  c(p)  , 

where  n is  the  minimal  number  of  steps  required  (so  that  h « n‘*  is  the  maximal  step- 
size  permitted),  and  the  cost  per  step  c(p)  is  the  number  of  arithmetic  operations 
required  for  the  execution  of  one  step  of  a p*h-order  NRK  method.  As  in  Traub  and 
Wofniakowski  [76],  we  shall  express  the  cost  per  step  in  the  form 

(3.3)  c(p)  e(Wp(v))  * d(p)  . 

Here  fllp(v)  is  the  information  about  v required  to  perform  one  step  of  a p^-order 
NRK  method  *p,  and  we  write  e($p(v))  for  the  informational  cost  of  *p;  we  call  d(p) 
the  combinatory  cost  of  ep.  For  example,  Euler’s  method 

Xj+i  - Xj  + h v(xj) 

has  informational  cost 

(3.4)  e(v)  s-  cost  of  evaluating  v at  one  point. 

The  combinatory  cost  is  two  operations  (i.e.,  one  addition  and  one  multiplication). 

We  now  assume  that  the  solution  x of  (1.1)  is  analytic  on  I.  Thus  Cauchy’s 
Integral  Theorem  (Ahlfors  [66],  pg.  122)  shows  that  there  exists  an  M > 0 such  that 


jx<k)(t)|  / K!  i M*  for  all  t « I . 

Finally,  we  shall  restrict  our  attention  to  problems  which  are  "sufficiently  difficult,"  i.e., 
for  which  there  exists  an  > 0 independent  of  h and  p so  that 

(3.5)  »G<Pp.h)  2 (Ml^P  if  h«  I and  p«  . 

(See  Section  4 of  Werschulz  [76b}) 

We  will  now  derive  a lower  bound  for  the  complexity  C(p,e)  via  NRK  methods. 
Clearly,  Theorem  2.1  implies  that  for  any  pth-order  NRK  method,  we  must  have 

(3.6)  e(Jlp(v))  2 e(v)  (p  + 1)  / 2 , 
and  a linear  lower  bound  on  the  combinatory  cost  states  that 

(3.7)  d(p)  2 aLp 

for  some  aL  > 0.  By  (3.6)  and  (3.7),  a lower  bound  on  the  cost  per  step  for  *p  is 

(3.8)  cL(p)  - (aL+e(v)/2)p  + e(v)/2 , 
which  leads  to 

Theorem  3.1: 

C(p,a)  2 Cj_(p,«)  :■  Ml  [(«l  + e(v)/2)  p + e(v)/2]  ea/P  . 

Proof:  From  (3.5),  we  see  that  if  (3.1)  holds,  then 

h i hL(p,«)  :■  Mf1  e"®/p  , 

Using  this  result,  (3.2),  and  (3.8),  the  theorem  follows.  | 

Next,  we  consider  upper  bounds  on  the  number  of  operations  required.  Instead 
of  using  ^ MSRK*  we  use  c*ass  *BRK  "unm°dified"  BRK  methods  described  in 
the  Appendix,  where  it  is  shown  that  4gp|(  is  order-convergent  In  the  sense  of 
Werschulz  [76b J.  That  is,  there  is  an  My  > 0 such  that 

(3.9)  *G(  Vh)  S (MU  h)P  » 

no  such  bound  is  Known  for  #mBRk>  In  addition,  ♦mbRK  r®duires  the  solution  of  p - 1 


? 


linear  systems  of  equations,  the  i*b  having  p - i unknowns,  in  order  to  perform  a 
“reorthogonalization."  So  the  smallest  known  combinatory  cost  for  this  class  is  about 
0(p3-81)  arithmetic  operations;  this  is  obtained  by  using  Strassen’s  technique  for  linear 
systems  (described  in  Borodin  and  Munro  [75]).  On  the  other  hand,  most  of  the 
combinatory  cost  for  pp  in  is  involved  in  finding  the  coefficients  of  the 

polynomial  pn+j  (see  the  Appendix);  once  these  coefficients  are  known,  the  remaining 
combinatory  cost  is  0(p  In  p)  as  p T oo.  An  estimate  of  how  much  work  is  required  to 
compute  these  coefficients  is  given  in 

Lemma  Let  xp,  yj, ... , yr,  w0,  Zp,  - • 2r  be  *iven>  and  ,et 

«»>  - *"o‘  «i 

be  the  unique  polynomial  of  degree  at  most  r 1 satisfying 

Q(xp)  “ w0  • Q'<*0>  " z0  » and  Q'ty)  " zi  <l  * j s r>  • 

If  T(r)  is  the  time  required  to  compute  qp, ... , qr+i,  then 

T(r)  - (Xr  ln2r)  as  r t oo  . 

Proof:  The  coefficients  qj,  •••  » (r-»-l)qr+j  of  Q'  may  be  computed  In  time 
0(r  ln2r)  by  using  a fast  algorithm  for  computing  the  coefficients  of  the  Lagrange 
polynomial  interpolating  the  points  (xp,  Zp),  (y^),  ...  , (yr,zrh  see  Borodin  and 
Munro  [75]  for  details.  Then  0(r)  operations  yield  qj, ... , qr+i,  and  Horner’s  rule  gives 
qp  with  0(r)  additional  operations.  | 

Thus  there  exists  a^  > 0 such  that 

(3.10)  d(p)  S «(jp  ln2(p+e)  . 

(We  write  "In  (p+e)",  where  e is  the  base  of  the  natural  logarithms,  rather  than  "In  p" 
as  a technical  convenience.  However,  an  expression  of  the  form  "In  (p+if)"  with  y > 0 
is  necessary  to  guarantee  that  d(l)  > 0.)  In  order  to  simplify  matters  a bit,  note  that 
Theorem  A.1  of  the  Appendix  implies  that 
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(3.11)  e(jRp(v))  £ e(v)  p . 

Although  tho  estimate  above  is  not  exact  for  p > 2,  it  is  asymptotically  equal  to  that  in 
Theorem  A.l.  (If  necessary,  the  sharper  estimate  given  there  may  be  used,  but  the 
calculation  of  optimal  order  (see  below)  involves  considerably  more  detail;  moreover, 
the  asymptotic  formulae  for  optimal  complexity,  order,  and  step-size  are  the  same  in 
either  case.)  Combining  (3.10)  and  (3.1 1),  we  see  that  the  cost  per  step  is  bounded  by 

(3.12)  cy(p)  - e(v)  p + ay  p ln2(p+e)  , 
which  leads  to 

Theorem  3.2: 

C(p,«)  £ Cy(p,a)  My  [e(v)  p + ay  p ln2(p+e)]  e“/p  . 

Proof:  If  we  set 

h - hy(p,«)  :■  My-1  e"°/p  , 

we  find  that  (3.9)  implies  that  (3.1)  holds.  Using  this  result,  (3.2),  and  (3.12),  the 
theorem  follows.  | 

Thus  we  have  found  bounds 

(3.13)  CL(p,«)  £ C(p,e)  £ Cy(p,«) 

on  the  number  of  operations  required  for  a p^-order  NRK  method  to  provide  an 
approximate  solution  satisfying  (3.1).  We  would  lihe  to  compute 

(3.14)  C*(«)  :-  inf  {C(p,o>  :p(  2 ++)  . 

This  is  not  possible,  since  we  only  have  bounds  for  C(p,e),  and  hence  cannot  compute 
C(p,«)  exactly.  However,  we  can  pick  optimal  orders  which  minimize  these  bounds. 
First,  we  prove 

Lemma  3.2:  Define 

GL(p)  s-  P2  Cl'(p)  / cL(p)  and  Gy(p)  :-  p2  cy'(p)  / Cy(p)  . 

Then  for  p > 0,  we  have  G|/(p)  > 0 and  Gy'fp)  > 0. 
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Proof:  Since  cl  is  a linear  polynomial  with  a negative  zero,  the  first  part  follows 
immediately.  Now  write  Cy(p)  - cj(p)  C2<p),  where 

cj(p)  :»  p and  C2<p)  :■  1 ♦ fi  ln2(p  ♦ e) , 
with  0 ay  / e(v) . Define 

Gj(p)  :-  p2  Cj'(p)  / c,(p)  (i  - 1,  2)  . 

Clearly  Gj'(p)  > 0 if  p > 0.  Now 

G2<p)  ■ 2 0 p2  In  (p+e)  / D2(p) , where  D2<p)  :■  (p+e)  • 

so  that 

G2/(p)  - 2 0 p g2(p)  / D2<p)2  , 

where 

g2<p)  :■  0 p ln2(p+e)  [In  (p+e)  - 1]  + 2 0 e ln2(p+e)  ♦ (p  + 2e)  In  (p+e)  ♦ p . 
Thus  G2/(p)  > 0 for  p > 0.  Since  Gy  - Gj  + G2  , the  desired  result  follows.  | 

We  now  have  the  following 

Theorem  3.3:  For  any  a > 0,  there  exist  P(_*(«)  and  py*(«)  such  that 
a - GL(p)  iff  p - Pl*<«)  and  • ■ G|j(p)  W P “ Pu*(«)  * 
Moreover, 

CL*(«)  :-  CL(pL*(«>,«)  < CL(p,«)  unless  p - pL%) 

and 

Cy*(«)  :-  Cyfpu^o),*)  < Cy(p,a)  unless  p - py*(«)  . 

Proof:  Using  (3.5),  (3.9),  and  Lemma  3.2,  this  follows  immediately  from  Lemma  2.1 
of  Werschulz  [76a].  | 

From  (3.13),  (3.14),  and  the  above  Theorem,  we  have  bounds 
(3.15)  CL*<«)  S C*<«)  S Cy*(«)  . 

We  call  pL*(a)  (respectively,  py*(«))  the  lower  (upper)  optimal  QUifiL  Cl*(«) 
(respectively,  Cy*(«))  the  lower  (upper)  optimal  complexity,  and 
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(3.16)  hL*(«)  hL(p|_*(«)^r)  (respectively,  hy^#)  hy(pu*(«),«)) 

the  lower  (upper)  optimal  step-size.  We  now  examine  how  these  quantities  behave  as 
a increases. 

Theorem  Pl*<«>.  Pu*<«>»  cL*<b)*  and  cU*(-,)  al1  iftCraa8a  monotonicelly  end 
tend  to  infinity  with  a.  Moreover,  the  following  asymptotic  formulae  hold  as  « tends  to 
infinity. 

(1.)  pj_*(«)  ~ a and  py*(«)  ~ « • 

(2.)  CL*(e)  ~ ML  e [aL  ♦ e(v)/2]  a and  Cy(«)  ~ My  ay  e « In2#  . 

(3.)  hL*(«)  -v  (Ml  e)'1  and  hy*<«)  ~ (Mye)-1  . 

Proof:  The  first  statement  follows  from  Lemma  3.2  and  from  Theorem  2.3  of 
Werschulz  [76b}  Now  Lemma  3.2  implies  that 

GL(p)  * P and  Gy(p)  ~ P as  p t oo  . 

Using  this  result  and  the  fact  that  lim  Pl*(«)  ■ liw  «too  * +00» 

(1.)  follows.  Finally,  (2.)  and  (3.)  follow  from  (1.),  Theorem  3.1,  and  Theorem  3.2.  | 

So  in  the  class  of  nonlinear  Runge-Kutta  methods,  we  find  that 

(3.17)  CL*(a)  - 0(«)  i C *(«)  S Cy*(«)  - 0(«  In2#) 
as  e tends  to  infinity;  so,  the  ratio 

Cy*(«)  / CL*<«)  - 0(ln2«)  as  • t oo 

indicates  the  gap  in  our  Knowledge  of  the  complexity  of  nonlinear  Runge-Kutta 
methods. 

Finally,  we  wish  to  compare  the  complexities  of  NRK  methods,  Taylor  series 
methods,  and  linear  Runge-Kutta  CLRK")  methods.  We  write  CytyRK*,  Cu,LRK*'  ^U,T* 
for  Cy*  in  the  class  of  NRK  methods,  LRK  methods,  and  Taylor  series  methods;  other 
notations  (CLLRK*,  CLRK*,  etc.)  are  formed  in  an  analogous  manner.  Finaly,  If 
f,  g : -»  satisfy  lim  .too  f(«)  - lim  bToo  g(«)  - ♦«,  we  write 
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<3.18)  f < g iff  Ha)  - o(g(«»  as  « t oo  | 

we  say  f is  asymptotically  less  than  g.  (See  Section  5 of  Werschulz  [76a].)  We  then 
have 

Theorem  3.5: 

CU,NRK*  < CL,LRK*  • 

(2.)  Cyj*  < Cy^RK  if  the  cost  of  evaluating  the  kth  derivative  of  v is 
bounded  for  all  k. 

Proof:  Immediate  from  (3.20)  and  (4.14)  of  Werschulz  [76a]  and  (3.17).  | 

As  a corollary  we  see  that  < Cj_pk*,  s0  that  the  best  NRK  method  known 

is  better  than  the  best  LRK  method  possible.  Moreover,  if  the  derivatives  of  v are 
easy  to  evaluate,  the  best  Taylor  series  method  known  is  better  than  the  best  NRK 
method  known.  However,  if  the  cost  of  evaluating  the  k*h  derivative  of  v increases 
faster  than  0(ln  k)  as  k t oo,  then  it  is  easy  to  show  thst  the  opposite  will  be  true. 


Appendix:  Order-Convergence  of  a Basic  Sequence 


In  this  Appendix,  we  describe  a subclass  of  a class  of  iterative  methods  for  tho 
solution  of  scalar  nonlinear  equations.  This  subclass  will  then  be  used  to  generate  an 
order-convergent  basic  sequence  tggg  of  nonlinear  Runge-Kutta  methods. 

Lemma  A.1:  Let  F:  DcR  -»  R have  a simple  zero  f,  and  suppose  that  F is 
analytic  at  f.  Pick  k,  m < Z+*  with  m ♦ 1 i k.  Then  there  is  a sequence 
^km  f^kmn  : n * Z -M‘}  of  stationary  multipoint  methods  without  memory  such  that 
the  following  hold: 

(1.)  The  method  ^mn  uses  the  information 

Wkmn(F)  {F(*0>'  "*F<mW»F(K)<yi>*  ”*F<K>(yn>> 

(the  points  yj  , ...  , yn  being  suitably  chosen)  to  compute  a new 
approximation  X|  to  f from  a given  approximation  xq  by  setting 

*1  *kmn<*0>  • 

(2.)  There  exists  a B > 0 and  an  hQ  > 0 such  that  if  |xq  - f | S hQ  , then 
|xj  - r|  £ (B  |x0  - f|)F  for  all  n< 

where 


(A.l)  p t"  min  (m  ♦ 2n  ♦ 1,  2m  ♦ n + 1)  . 

Before  proving  the  Lemma,  we  describe  how  the  method  computes 

improved  approximation  xj  from  the  old  approximation  xq  . 
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Algorithm  for  computing  xt  ^rW 
(1.)  Lot  8 |F(xq)/F'(xq)| . 

(2.)  Lot  bo  an  approximate  zero  of 

Pl<*>  <x  " x0>‘  / i! 

aatisfying 

(A.2)  zi  - xq  ♦ 0(0  and  |p1(z1>|  S (AA  l)m+1  , 
where  Aj  is  independent  of  n. 

(3.)  Let 

V\  im  *0  * *in  <Z1  ' X0J  <l  * 1 * n> » 

where 

ain  ^ * xin^  I ^ 

and  xln  > ...  > xnn  are  the  zeros  of  the  Jacobi  polynomial 
Pn(x)  Pff-1' m+l'K)(x) 

(see  SzegB  [59]). 

(4.)  Let  pn+j  be  the  polynomial  of  degree  at  most  m ♦ n that  interpolates  the 
information  Whmn^»  arK<  *•*  X1  #n  •PProx'm**a  z#r0  °*  Pn+1 
satisfying 

(A.3>  Xi  - x0+CXO  and  |pn+i(xi>l  i (A2  •>* . 
where  A2  is  independent  of  n and  p is  given  by  (A.1). 

Here  we  use  the  notation  of  Brent  [74].  Clearly,  < C#(k,  m,  n),  the  only 
difference  being  that  conditions  (A.2)  and  (A.3)  replace  (2.2)  and  (2.4)  of  Brent  [7 4}.  It 
is  easy  to  see  that  (A.2)  and  (A.3)  may  be  realized  by  using  riog2(m+l)1  - l and 
riog2(e/<m+l))1  iterations  of  Newton’s  method,  with  the  respective  starting 
approximations  of  xq  - F(xq)  / F'(xq)  and  zj  . 


Proof  of  Lemma  A.1:  Let  x^'  be  the  exact  zero  of  pn+j  near  xq.  We  then  find 
that  there  is  a ( between  xj'  and  zj  such  that 

(A.4)  IRxjOl  i IPn*l(*l>  ’ F(*l>l  ♦ W+l<t>  * r'<M  IV  *ll  • 

Using  (A.3),  the  analyticity  of  F,  and  standard  techniques  of  interpolation  theory 
(Traub  [64]),  it  is  easy  to  show  that  (2.9)  and  (2.10)  of  Brent  [74]  may  be  rewritten  as 

|pn+l(x)  - F(x)|  S (A3  l)m*"+l  and 

(A.5) 

|p^+l(x)  - F'(x)|  i (A4l)m+n 

for  |x  - xq|  s 41.  (Here  all  constants  Ar  will  be  independent  of  n.)  Similarly,  we  find 
that 

' - ri  * <a5  •>"**"  «*■  i*i  • fl  * (a6  ,>m*1  • 

so  that  the  triangle  inequality  gives 

(A.6)  |x!' -Zj|  S (A7l)m+1  . 

Using  (A.4),  (A.5),  and  (A.6),  we  see  that 

|f<xl0|  S |pn+l(*i>  - Rzj)|  ♦ <A8  |)2m+n+i 

S |pn+l<*!)  - F^H  ♦ |F2(zj)|  ♦ (Ag  , 


where 


F^x)  X-2"  (x  - x0)'  F<')(x0)  / i!  and  F2(x)  F(x)  - F^x)  . 


Clearly  |F2(x)|  S (Ag  |)m+2n*l  ( 80  that  (A.7)  becomes 

(A.8)  |F(Xl')|  S lPn+l^zl>  “ Fi(*i>l  ♦ (A10  *)F  • 

As  in  Brent  [74],  we  now  write 

pn+i<x)  ■ r^(x)  ♦ r2(x) , 

where  rj  (i  - 1,  2)  is  the  polynomial  of  degree  at  most  m ♦ n satisfying 

rjW(XQ>  - F}>X*0)  (0  & j S m) 


rf(K)(yj)  - Fj^(yj)  (l  S j i n)  . 
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If  we  let 


Pix)  r j(x  ♦ Xq)  - Fj(x  ♦ xq)  , 
end  write  a zj  - xq  (in  this  Appendix  only),  we  find  that 

p<')(0)  - 0 (0  S i 5 m)  and  P<K)(ejn  •)  - 0 (1  S I i n)  . 
We  may  easily  alter  the  proof  of  Lemma  4.3  in  Brent  [74]  to  show  that 

r !<*!>-  F^)  - P(i)  - 0 . 

Thus  (A.8)  becomes 


(A.9)  |F(xj')|  S ♦ (A10  S)A. 

To  bound  the  remaining  term,  let  us  write 

r2(x>  - XjLj  sj+(n  (x  - x0)i+m  , 

recalling  that  r2  has  a zero  of  multiplicity  m at  xq  . Using  the  notation  of  Stewart  [73], 
we  see  that  the  nonzero  coefficients  of  are  given  by  the  solution  of  the  linear 
system 


Wy  - c , 

where 

«,j  «|’1  (1  S i,  j & n) , 

aj+m  (j  ♦ m)»  / (j  ♦ m - k)l  (1  i j S n) , and 
T,  «K  F2<K><yi)  / «™'K+1  (1  S i S n)  . 

Since  WT  is  a Vandermonde  matrix,  we  find  that  the  entries  of  U • W“*  are  given  by 
tfij  " *jn  ' *n-i,n-l,j  I ^*jn  ” *rn^  * 

where 


V"-l»j  •"*  1 *Pi»n  "•  *Pj,,n  • 

the  sum  being  taken  over  all  multi-indices  Pf  ...  p^  not  including  j (Gregory  and 
Karney  [69]).  Since  there  are  fewer  than  2n  summands,  each  of  which  lies  in  [0,  1] , 
we  see  that  * 2n  , implying  that 


(A.10) 


(see  Abramowitz  and  Stegun  [64]). 


Now  it  is  clear  that 

, . m-k 

lSjSn  1 / «jn  - 1 
By  Theorem  8.9.1  of  Szegfi  [59],  we  may  show  that 


using  this  result  and  (22.5.2)  of  Abramowitz  and  Stegun  [64],  we  find  that 

lsjsn  t«£’K  Gn/<«jn>3"1 

<A.l  1>  + 

i A12  n2<m-R)  ( 1 max  lsjsn  |Pn'<Xjn>r 


By  the  symmetry  relation  (4.1.3)  of  SzegO  [59],  we  may  assume  that  0 i x,-n  < 1.  Using 
Theorem  8.9.1  of  Szeg6  [59],  we  may  show  that 

ip„'(,jnr‘  s i3)" . 

and  so  (A.10),  (A.  11),  the  definition  of  F2,  and  the  above  imply  that 

1.1  < *a.  . i\m+2n*l 


yielding  the  result 


,m+2n+l 


So  (A.9)  becomes 


By  Taylor's  Theorem,  this  implies 


19 


I*!'  - f|  i (Aj7  W . 

The  desired  result  then  follows  from  (A.3)  and  from  (2.5)  of  Brent  [74].  | 

We  now  describe  the  basic  sequence  tggg  . The  methods  in  this  basic  sequence 
are  given  by 

di<XQ.W  *-  v(k0), 

*2<*0  » h>  !"  ^x0  + h v(x0*  f • 

and  for  p 2 2, 

Pp(*0  * h)  !"  h_l  1*1,1, p-2(x0)  ’ x0l  • 

with  *°  the  function  F given  by  (2.8)  and  the  approximation  xj  to  Xj' 

being  given  by  an  appropriate  number  of  iterations  of  Newton’s  method  (as  described 
above). 

Theorem  A.1:  The  basic  sequence  tgpiC  '*  order-convergent  with  respect  to  the 
global  error.  Moreover,  the  number  of  stages  s(p)  required  by  pp  « 9gRK  '*  **ven  by 


Proof:  We  use  the  notation  of  Lemma  A.1,  writing  z(h)  for  the  computed  p^- 
order  approximation  xj  to  x(h)  and  pn+j(  ■ , xq)  for  the  polynomial  pn+j  . The  result 
of  Lemma  A.1  is  that 

h’1  |z(h)  - x(h)|  S (B  h)P  , 

the  desired  result  for  a single  unit  step.  To  prove  the  global  result,  we  must  consider 
the  Lipschitz  constants  for  tgRK* 

We  implicitly  differentiate  the  result  Pn+^xj'i  xq)  ■ 0 to  find 

d\  pp(*o  • h>  " "h-1  <Wxl''  x0>  * ,p<x0> ' 

where 

<W*l'.  x0>  “ 1 ♦ d2  Pn*l<xl'»  x0>  / dl  Prt+l<xl'»  x0> 
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and 

•p(x0)  - h"1  (d/dx0)  [xj  - Xj']  . 

It  is  easy  to  see  that  xj  and  xj'  are  analytic  functions  of  xg  . Since  their  difference 
tends  to  zero  uniformly  on  the  domain  of  v as  p t oo  , it  follows  that 

,im  ptco  «p<*0>  " 0 • 

We  claim  that 

Qn+jfxi',  x0)  - 0(h  In  n)  as  n t oo  , 

uniformly  in  xg  . To  see  this,  note  that  we  may  vrite  the  interpolation  polynomial  pn+i 
in  terms  of  Jacobi  polynomial  Pn  , finding  that 

Pn*l<*.*0>  " (-0°  <h/2>  $[{*}  Pn(t)  dt  ♦ h v(x0)  Ihn  - h, 

where 

f(x)  2 (x  - Xg)  / [h  v(Xg)]  - 1 

and 

■kn  (2  « * W v<yk)  Pn'O^r*  « * » P„«>  / « - W <*  • 

Now 

- <-l)nPn<ri)/v(x0)  ♦ (1  ♦ f 1)  I «<*kn>  lWfl> . 

where 

r l 

!"  ^n^  f P*n^xkn)  ^x  " xkn^  » ,n<* 
g(t)  1 / t(l  ♦ t)  v(x0  ♦ (1  ♦ t)  h vfxg)  / 2)]  . 

By  (8.2 1. 10)  of  Szegd  [59],  the  first  term  in  the  expression  for  dj  Pn+l(*l' » x0^  I0** 
to  zero  as  n T oo  . A minor  modification  of  the  proof  of  Theorem  14.4  of  SzegO  [59] 
shows  that  the  sum  in  the  remaining  term  tends  to  g(f(x(h)))  as  n t oo  . So 

dl  Pn+l*xl'  • * v<x<h))"1  «*  n * 00  • 


' 


Using  Lemma  A.1  of  the  Appendix  in  Werschulz  [76a]  and  techniques  similar  to  those 
yielding  the  above  estimate,  we  find 

^2  Pn+l^l'  * *(P  " 0(h  In  n)  - v(x(h))"1  as  n t co  . 

This  gives  the  estimate  claimed  for  Qn+i(xi' , x0) . 

So  the  Lipschitz  constant  for  < %gRK  6rOWS  ••  logarithm  of  p.  By 
Proposition  4.3  of  Werschulz  [76b],  4ggg  is  order-convergent.  | 


i; 
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